Plantear el problema de minimos cuadrados como un problema de minimizacion:
Tenemos que \(Y=X\beta + \epsilon\), por lo que donde \(\epsilon\) es el error, por lo que podemos plantear el problema como la minimizaci´ón de este error \[min_{ } ||\epsilon||^2 \] que es equivalente a \[min_{\beta} ||Y-X\beta||^2\] Para encontrar el estimador \(\hat{\beta}\) es necesario derivar e igualar a 0: Derivando respecto a \(\beta\): \[-2X^t (Y-X\hat{\beta})=0\] \[2X^tY=2X^tX\hat{\beta}\] \[\hat{\beta}=(X^tX)^{-1}X^tY\] Es el \(\hat{\beta}\) que soluciona el problema de minimización.
Con esto podemos ver que para cada: \[\hat{\beta}_{i}=\frac{\sum_{i}x_{i}y_{i}}{\sum_{i}x_{i}^2}\] La cuál es una composición lineal de los parametros. Desde el planteamiento estamos buscando que el resultado sea lineal en los parametros, no en los datos por lo que este metodo serviria para poder aproximar polinomios de la orma \(Y=X^2\)
Recordemos que al encontrar la proyección de un vector \(Proy_{X}(y)= \frac{<x,y>}{||x||^2}y\) el error de proyeccion \(Y- Proy_{x}(Y)\) será ortogonal a la proyección. En el caso del problema de minimizar el error, nosotros estamos proyectando la variable dependiente Y sobre el espacio formado por nuestros datos X. y buscamos las combinaciones sobre X que hagan este error mas pequeño.
Viendolo como el teorema de Pitagoras el cu´ál nos dice que el cuadrado de la hipotenisa de un triangulo es igual a la suma de los otros dos lados \(h^2=x^2+y^2\) si el angulo formado entre \(x\) y \(y\) es de 90 grados, en este caso \(H=Y\) las variables dependientes, \(x=Proy_{x}(Y)\) y \(y= error\), entonces la manera en la que se puede cumplir el teorema es que el error sea ortogonal a la proyeccion (angulo de 90 grados) sino el error seráa mas grande.
Al definir una columna de unos estamos cambiando un poco el problema. Antes de agregar la columna: \[y_{j}=\beta_{1}x_{1j}+...+\beta_{n}x_{nj}\] Al definir la columa de unos estamos dando cierta holgura al modelo ya que al agregar un parametro \(\alpha\) no estamos forzando a que la aproximacion pase por el origen. \[y_{j}=\alpha+\beta_{1}x_{1j}+...+\beta_{n}x_{nj}\]
Planteando el problema como \[Y= X\beta +\epsilon\] con \[\epsilon \sim N(0,\sigma^2I_{n})\]
podemos probar que \[Y \sim N(X\beta,\sigma^2I_{n})\]
Demostracion: \[Y= X\beta +\epsilon\] \[\epsilon=Y-X\beta\]
Ahora, falta demostrar que Y es Normal: Por propiedades de la Normal como \(\epsilon \sim N(0,\sigma^2I_{n})\) \[\epsilon + X\beta \sim N(X\beta, \sigma^2I_{n})\] (Esto era m´ás rapido pero bueno…)
Entonces con esto podemos escribir la funcion de verosimilitud de Y como: \[ L(\beta, \sigma^2)=\prod \frac{1}{\sqrt{2\pi\sigma^2}} exp(-\frac{(y-X\beta)^2}{2\sigma^2})\]
\[=(\frac {1}{\sqrt{2\pi\sigma^2}})^n exp (- \frac{1}{2\sigma} (y-X\beta)^t(y-X\beta))\]
A la cual le sacamos Logaritmo para despu´és poder buscar los parametros que maximicen la funcion: \[Log(L(\beta, \sigma^2))=-\frac{n}{2}ln(2\pi)-\frac{n}{2}ln(\sigma^2)-\frac{1}{2\sigma^2}(y-X\beta)^t(y-X\beta)\] Ahora hay que maximizar esa funcion sobre \(\beta\) y \(\sigma^2\) Respecto de \(\beta\) \[\frac{dL}{d\beta}=\frac{1}{\sigma^2}(y-X\beta)^tX=0\] \[X^t y=X^tX\beta\] \[\hat{\beta}=(X^tX)^{-1}X^ty\]
Respecto de \(\sigma^2\) \[\frac{dL}{d\sigma^2}= - \frac{n}{2\sigma^2}+\frac{1}{2\sigma^4}(y-X\beta)^t(y-X\beta)=0\] \[\sigma^2=\frac{1}{n} (y-X\beta)^t(y-X\beta)\] que puede resolverse usando \(\hat{\beta}\) que encontramos en la solucion anterior.
Como pudimos ver el resultado de maxima verosimilitud \[\hat{\beta}=(X^tX)^{-1}X^ty\] Es igual al resultado de \(\hat{\beta}\) de minimos cuadrados
El teorema de Gauss Markov plantea que un modelo de regresion lineal en el que se cumple que:
Entonces el mejor estimador lineal insengado de los oeficientes \(\beta\) es el resultante del problema de minimizar los errores al cuadrado. si existe.
library(ggplot2)
library(plotly)
data("diamonds")
head(diamonds)
#Entonces, para hacer la regresion lineal tengo que agarrar los #valores numericos que son
diamonds_num<-diamonds[,c(1,5:10)]
reg_precio=lm(price~carat+depth+table+x+y+z, data=diamonds_num)
summary(reg_precio)
Call:
lm(formula = price ~ carat + depth + table + x + y + z, data = diamonds_num)
Residuals:
Min 1Q Median 3Q Max
-23878.2 -615.0 -50.7 347.9 12759.2
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 20849.316 447.562 46.584 < 2e-16 ***
carat 10686.309 63.201 169.085 < 2e-16 ***
depth -203.154 5.504 -36.910 < 2e-16 ***
table -102.446 3.084 -33.216 < 2e-16 ***
x -1315.668 43.070 -30.547 < 2e-16 ***
y 66.322 25.523 2.599 0.00937 **
z 41.628 44.305 0.940 0.34744
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1497 on 53933 degrees of freedom
Multiple R-squared: 0.8592, Adjusted R-squared: 0.8592
F-statistic: 5.486e+04 on 6 and 53933 DF, p-value: < 2.2e-16
Que tan bien ajusta el modelo depende de la correlacion que tenga cada una de las variables explicativas sobre el precio, esto lo podemos ver a trav´és de:
diamonds_num<- diamonds_num[,c(4,1,2,3,5,6,7)]
cor(diamonds_num)
price carat depth table x
price 1.0000000 0.92159130 -0.01064740 0.1271339 0.88443516
carat 0.9215913 1.00000000 0.02822431 0.1816175 0.97509423
depth -0.0106474 0.02822431 1.00000000 -0.2957785 -0.02528925
table 0.1271339 0.18161755 -0.29577852 1.0000000 0.19534428
x 0.8844352 0.97509423 -0.02528925 0.1953443 1.00000000
y 0.8654209 0.95172220 -0.02934067 0.1837601 0.97470148
z 0.8612494 0.95338738 0.09492388 0.1509287 0.97077180
y z
price 0.86542090 0.86124944
carat 0.95172220 0.95338738
depth -0.02934067 0.09492388
table 0.18376015 0.15092869
x 0.97470148 0.97077180
y 1.00000000 0.95200572
z 0.95200572 1.00000000
Como podemos ver en la tabla de correlaciones, precio est´á correlacionada altamente con “carat”, “x”,“y”,“z” y tienen cierta correlacion con death y table, esto me dice que las variables escogidas sirven para explicar el precio. Podemos ver esta relacion de manera grafica a continuacion en el siguiente Diagrama de dispersion :
library(GGally)
package ‘GGally’ was built under R version 3.2.5replacing previous import by ‘utils::capture.output’ when loading ‘GGally’replacing previous import by ‘utils::head’ when loading ‘GGally’replacing previous import by ‘utils::installed.packages’ when loading ‘GGally’replacing previous import by ‘utils::str’ when loading ‘GGally’
library(ggplot2)
my_fn <- function(data, mapping, ...){
p <- ggplot(data = data, mapping = mapping) +
geom_point() +
geom_smooth(method=loess, fill="red", color="red", ...) +
geom_smooth(method=lm, fill="blue", color="blue", ...)
p
}
g = ggpairs(diamonds_num,columns = 1:7, lower = list(continuous = my_fn))
g
plot: [1,1] [=--------------------------------] 2% est: 0s
plot: [1,2] [=--------------------------------] 4% est: 5s
plot: [1,3] [==-------------------------------] 6% est: 5s
plot: [1,4] [===------------------------------] 8% est: 5s
plot: [1,5] [===------------------------------] 10% est: 5s
plot: [1,6] [====-----------------------------] 12% est: 5s
plot: [1,7] [=====----------------------------] 14% est: 4s
plot: [2,1] [=====----------------------------] 16% est: 4s
plot: [2,2] [======---------------------------] 18% est: 2m
plot: [2,3] [=======--------------------------] 20% est: 2m
plot: [2,4] [=======--------------------------] 22% est: 2m
plot: [2,5] [========-------------------------] 24% est: 1m
plot: [2,6] [=========------------------------] 27% est: 1m
plot: [2,7] [=========------------------------] 29% est: 1m
plot: [3,1] [==========-----------------------] 31% est: 1m
plot: [3,2] [===========----------------------] 33% est: 2m
plot: [3,3] [===========----------------------] 35% est: 3m
plot: [3,4] [============---------------------] 37% est: 2m
plot: [3,5] [=============--------------------] 39% est: 2m
plot: [3,6] [=============--------------------] 41% est: 2m
plot: [3,7] [==============-------------------] 43% est: 2m
plot: [4,1] [===============------------------] 45% est: 2m
plot: [4,2] [===============------------------] 47% est: 2m
plot: [4,3] [================-----------------] 49% est: 2m
plot: [4,4] [=================----------------] 51% est: 3m
plot: [4,5] [==================---------------] 53% est: 2m
plot: [4,6] [==================---------------] 55% est: 2m
plot: [4,7] [===================--------------] 57% est: 2m
plot: [5,1] [====================-------------] 59% est: 2m
plot: [5,2] [====================-------------] 61% est: 2m
plot: [5,3] [=====================------------] 63% est: 2m
plot: [5,4] [======================-----------] 65% est: 2m
plot: [5,5] [======================-----------] 67% est: 2m
plot: [5,6] [=======================----------] 69% est: 2m
plot: [5,7] [========================---------] 71% est: 2m
plot: [6,1] [========================---------] 73% est: 2m
plot: [6,2] [=========================--------] 76% est: 2m
plot: [6,3] [==========================-------] 78% est: 2m
plot: [6,4] [==========================-------] 80% est: 1m
plot: [6,5] [===========================------] 82% est: 1m
plot: [6,6] [============================-----] 84% est: 1m
plot: [6,7] [============================-----] 86% est: 1m
plot: [7,1] [=============================----] 88% est: 1m
plot: [7,2] [==============================---] 90% est:48s
plot: [7,3] [==============================---] 92% est:40s
plot: [7,4] [===============================--] 94% est:31s
plot: [7,5] [================================-] 96% est:22s
plot: [7,6] [================================-] 98% est:11s
plot: [7,7] [=================================]100% est: 0s
Graficamente así se ven los valores del modelo contra los observados. Podemos ver que no forman una linea de 45 grados por lo que posiblemente hay alg´ún problema dada la correlacion entre las variables explicativas. Esto tambien lo podemos ver en la tabla de correlaciones.
#Grafica de los valores del modelo contra los reales
library(plotly)
data<-as.data.frame(cbind(reg_precio$fitted.values,diamonds_num$price))
p <- plot_ly(data = data, x = ~data[,1], y = ~data[,2]) %>%
layout(title="Valores ajustadosvs Valores observados")
replacing previous import by ‘shiny::includeHTML’ when loading ‘crosstalk’replacing previous import by ‘shiny::knit_print.shiny.tag’ when loading ‘crosstalk’replacing previous import by ‘shiny::code’ when loading ‘crosstalk’replacing previous import by ‘shiny::includeScript’ when loading ‘crosstalk’replacing previous import by ‘shiny::includeMarkdown’ when loading ‘crosstalk’replacing previous import by ‘shiny::tags’ when loading ‘crosstalk’replacing previous import by ‘shiny::is.singleton’ when loading ‘crosstalk’replacing previous import by ‘shiny::withTags’ when loading ‘crosstalk’replacing previous import by ‘shiny::img’ when loading ‘crosstalk’replacing previous import by ‘shiny::tagAppendAttributes’ when loading ‘crosstalk’replacing previous import by ‘shiny::knit_print.shiny.tag.list’ when loading ‘crosstalk’replacing previous import by ‘shiny::knit_print.html’ when loading ‘crosstalk’replacing previous import by ‘shiny::tagAppendChild’ when loading ‘crosstalk’replacing previous import by ‘shiny::includeCSS’ when loading ‘crosstalk’replacing previous import by ‘shiny::br’ when loading ‘crosstalk’replacing previous import by ‘shiny::singleton’ when loading ‘crosstalk’replacing previous import by ‘shiny::span’ when loading ‘crosstalk’replacing previous import by ‘shiny::a’ when loading ‘crosstalk’replacing previous import by ‘shiny::tagList’ when loading ‘crosstalk’replacing previous import by ‘shiny::strong’ when loading ‘crosstalk’replacing previous import by ‘shiny::tag’ when loading ‘crosstalk’replacing previous import by ‘shiny::p’ when loading ‘crosstalk’replacing previous import by ‘shiny::validateCssUnit’ when loading ‘crosstalk’replacing previous import by ‘shiny::HTML’ when loading ‘crosstalk’replacing previous import by ‘shiny::h1’ when loading ‘crosstalk’replacing previous import by ‘shiny::h2’ when loading ‘crosstalk’replacing previous import by ‘shiny::h3’ when loading ‘crosstalk’replacing previous import by ‘shiny::h4’ when loading ‘crosstalk’replacing previous import by ‘shiny::h5’ when loading ‘crosstalk’replacing previous import by ‘shiny::h6’ when loading ‘crosstalk’replacing previous import by ‘shiny::tagAppendChildren’ when loading ‘crosstalk’replacing previous import by ‘shiny::em’ when loading ‘crosstalk’replacing previous import by ‘shiny::div’ when loading ‘crosstalk’replacing previous import by ‘shiny::pre’ when loading ‘crosstalk’replacing previous import by ‘shiny::htmlTemplate’ when loading ‘crosstalk’replacing previous import by ‘shiny::suppressDependencies’ when loading ‘crosstalk’replacing previous import by ‘shiny::tagSetChildren’ when loading ‘crosstalk’replacing previous import by ‘shiny::includeText’ when loading ‘crosstalk’replacing previous import by ‘shiny::hr’ when loading ‘crosstalk’
p
No trace type specified:
Based on info supplied, a 'scatter' trace seems appropriate.
Read more about this trace type -> https://plot.ly/r/reference/#scatter
No scatter mode specifed:
Setting the mode to markers
Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
No trace type specified:
Based on info supplied, a 'scatter' trace seems appropriate.
Read more about this trace type -> https://plot.ly/r/reference/#scatter
No scatter mode specifed:
Setting the mode to markers
Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
Podemos ver que como el modelo tiene una \(R^2 ~ .8592\) por lo que podemos decir que el modelo explica 85% de la varianza con las variables que estamos usando pero hay algunos problemas entre ellas por lo que los estimadores deben de estar sesgados. nuestro error estandar residual \(\sigma=1497\) por lo que \(sigma^2=1497^2\)
\[arcos(R^2)=\theta\] En este caso en particular \(R^2=.8592\) por lo que
acos(sqrt(.8592))
[1] 0.3846484
En caso de que podamos aproximar con una normal tenemos lo siguiente:
lik1<-function(par,X,y){
beta<-par[1:ncol(X)]
sigma2<-par[ncol(X)+1]
n<-nrow(y)
X<-as.matrix(X)
beta<-as.matrix(beta)
mu<-X%*%beta
logl= -.5*n*log(2*pi) -.5*n*log(sigma2) - (sum((y-mu)**2)* .5/sigma2))
Error: unexpected ')' in:
"
logl= -.5*n*log(2*pi) -.5*n*log(sigma2) - (sum((y-mu)**2)* .5/sigma2))"
En nuestro caso en particular:
optim(par=betasigma,lik1, X=X_reg, y=y_reg,method="BFGS")
NaNs producedNaNs producedNaNs produced
$par
[1] 27260.5528 10906.0130 -271.6052 -134.3417 -1586.6727 271.2380 140486.4120
$value
[1] 802048.1
$counts
function gradient
159 100
$convergence
[1] 1
$message
NULL
reg_precio$coefficients
(Intercept) carat depth table x y z
20849.3164 10686.3091 -203.1541 -102.4457 -1315.6678 66.3216 41.6277
Ahora, como z es “no significativa” vamos a sacarla del modelo y a correr todo con el nuevo conjunto de variables
summary(reg_red)
Call:
lm(formula = price ~ carat + depth + table + x + y, data = diamonds_red)
Residuals:
Min 1Q Median 3Q Max
-23872.1 -614.8 -50.5 347.5 12759.4
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 20702.947 419.575 49.343 < 2e-16 ***
carat 10686.707 63.199 169.095 < 2e-16 ***
depth -200.718 4.855 -41.344 < 2e-16 ***
table -102.490 3.084 -33.234 < 2e-16 ***
x -1293.542 36.063 -35.869 < 2e-16 ***
y 69.575 25.287 2.751 0.00594 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1497 on 53934 degrees of freedom
Multiple R-squared: 0.8592, Adjusted R-squared: 0.8592
F-statistic: 6.583e+04 on 5 and 53934 DF, p-value: < 2.2e-16
Ahora con esto vuelvo a definir mis variables a usar en la funcion de max verosimilitud
optim(par=betasigma,lik1, X=X_reg, y=y_reg, method = "BFGS")
NaNs producedNaNs producedNaNs produced
$par
[1] 27260.5528 10906.0130 -271.6052 -134.3417 -1586.6727 271.2380 140486.4120
$value
[1] 802048.1
$counts
function gradient
159 100
$convergence
[1] 1
$message
NULL
Podemos ver como los estimadores obtenidos por maxima verosimilitud se parecen mas a los que tienen mayor significancia en el modelo de regresion lineal.